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Statement  of  Problem  and  Most  Important  Results 


Forward 

This  report  describes  the  research  accomplishments  of  the  project  “Scaling  Properties  and  Spatial 
Interpolation  of  Soil  Moisture”  conducted  at  Pennsylvania  State  University.  The  planned 
duration  of  this  project  was  five  years,  but  the  principal  investigator  relocated  from  Pennsylvania 
State  University  to  Colorado  State  University  after  approximately  one  year.  Consequently,  the 
grant  at  Pennsylvania  State  University  was  tenninated  in  May  2004  after  one  year  of  funding, 
and  a  new  grant  was  issued  at  Colorado  State  University  for  the  remaining  four  years  of  the 
project.  This  final  report  describes  the  progress  made  during  the  first  year  of  the  five-year 
project. 


1  Problem  Statement 

Soil  moisture  patterns  play  an  important  role  in  numerous  hydrologic  and  geomorphic 
applications.  They  directly  impact  many  hydrologic  processes  including  infiltration,  surface 
runoff,  transpiration,  and  recharge.  Through  these  processes,  soil  moisture  patterns  also 
influence  precipitation  recycling  within  a  region,  thus  perpetuating  wet  and  dry  periods  (Findell 
and  Eltahir,  1997).  Soil  moisture  also  impacts,  either  directly  or  indirectly,  geomorphic 
processes  such  as  bedrock  weathering  and  sediment  detachment,  transport,  and  deposition. 

Despite  these  important  roles,  soil  moisture  patterns  and  dynamics  are  not  well 
understood.  The  difficulty  in  understanding  soil  moisture  arises  from  two  sources.  First,  soil 
moisture  patterns  and  dynamics  are  very  complex,  and  variability  occurs  over  a  wide  range  of 
spatial  and  temporal  scales.  At  large  spatial  scales,  soil  moisture  depends  largely  on  climatic 
characteristics  and  soil  types,  but  at  smaller  scales,  topography  becomes  more  important  in 
determining  soil  moisture  patterns.  Furthermore,  the  roles  of  these  physical  characteristics  can 
vary  in  time.  During  dry  periods,  soil  type  is  more  important  in  determining  soil  moisture 


patterns,  whereas  in  wet  periods,  the  role  of  topography  is  enhanced  (Grayson  et  ah,  1997; 
Western  et  ah,  2001).  Soil  moisture  patterns  and  dynamics  also  depend  on  the  depth  over  which 
the  moisture  is  considered.  Near  the  surface,  soil  moisture  is  highly  dynamic,  but  it  becomes  less 
dynamic  with  increasing  depth.  The  proximity  of  the  water  table  also  impacts  the  nature  of  the 
variability  with  depth  (Salvucci  and  Entekhabi,  1994a,  1994b). 

The  second  difficulty  in  characterizing  soil  moisture  is  measurement  (see  Robock  et  ah, 
2000).  Numerous  in-situ  methods  are  available  to  characterize  soil  moisture  (neutron  probes,  for 
example),  but  such  methods  quantify  soil  moisture  only  at  the  sampling  locations.  Due  to  the 
large  and  complex  spatial  variability  observed,  such  measurements  are  not  easy  to  extrapolate  to 
the  spatial  scales  of  most  hydrologic  and  geomorphic  applications.  Alternatively,  remote  sensing 
methods  can  be  used  to  measure  soil  moisture,  either  by  aircraft  or  satellite.  Most  such  methods 
rely  on  passive  microwave  measurements.  Unfortunately,  remote  sensing  methods  characterize 
soil  moisture  only  near  the  ground  surface  and  provide  measurements  at  a  relatively  coarse 
resolution. 

2  Project  Objectives 

The  long-tenn  objective  of  this  research  project  is  to  develop  an  interpolation  or  downscaling 
method  for  soil  moisture  observations  that  accounts  for  the  role  of  topography.  To  accomplish 
this  primary  objective,  several  secondary  goals  must  be  achieved.  First,  a  better  understanding 
of  soil  moisture  patterns  is  required,  which  includes  characterizing  the  spatial  and  temporal 
variability  as  well  as  their  dependence  on  the  scale  of  observation.  Ultimately,  one  would  like 
the  ability  to  simulate  the  soil  moisture  variability  in  space  and  time  using  climatic  inputs. 
Second,  the  dependence  of  soil  moisture  patterns  on  a  region’s  physical  characteristics  must  be 
analyzed.  This  goal  involves  describing  how  soil  moisture  patterns  depend  on  variables  such  as 
soil  type  and  topography  at  different  scales.  Ultimately,  the  connections  between  topographic 
and  soil  moisture  scaling  characteristics  are  sought.  Third,  a  better  understanding  of  the  scaling 
characteristics  of  topography  and  their  dynamic  origins  is  needed.  In  particular,  one  would  like 
to  quantify  the  relationship  between  the  subsurface  redistribution  of  water  and  geomorphic 
processes  and  river  basin  form.  In  the  end,  this  understanding  will  guide  the  development  of  soil 
moisture  interpolation  procedures  that  include  the  role  of  topography. 


3  Approach 

The  general  approach  used  in  this  project  was  a  combination  of  statistical  analysis  of  digital  data, 
numerical  and  analytical  modeling  of  soil  moisture,  and  numerical  modeling  of  topographic 
dynamics.  Soil  moisture  data  from  three  sources  was  used:  the  Illinois  State  Water  Survey  data 
set,  the  Southern  Great  Plains  experiment  in  1997  (SGP  97),  and  the  Washita  experiment  in  1992 
(Robock  et  al.,  2000).  The  Illinois  data  was  collected  with  an  in-situ  method  (neutron  probes) 
and  describes  the  soil  moisture  with  only  crude  spatial  detail.  However,  this  dataset  has  a  long 
duration  (approximately  20  years),  so  it  is  well-suited  for  detailed  studies  of  soil  moisture 
dynamics.  The  SGP  97  data  and  the  Washita  data  were  primarily  collected  with  aircraft-based 
remote  sensing  methods,  but  in-situ  methods  were  used  to  ground-truth  the  data.  These  datasets 
have  good  spatial  resolution,  but  limited  duration.  The  experiments  collected  data  over  a  period 
of  a  few  days  to  weeks.  The  Washita  dataset  describes  a  subset  of  the  SGP  97  field  site.  All  of 
these  datasets  describe  regions  with  semi-humid  temperate  climates  and  low  relief. 

Three  primary  methods  of  statistical  analysis  were  used.  First,  basic  statistical  properties 
of  the  soil  moisture  patterns  were  measured.  This  approach  largely  neglects  spatial  and  temporal 
correlations  and  connectivity  for  simplicity.  Second,  Empirical  Orthogonal  Functions  (EOFs) 
were  extracted  from  the  soil  moisture  data  and  compared  with  ancillary  data.  This  approach 
allows  one  to  identify  the  dependence  of  soil  moisture  on  other  physical  characteristics  of  the 
region.  Third,  scaling  analyses  were  used  to  examine  how  the  basic  properties  of  soil  moisture 
change  with  the  scale  of  observation. 

The  modeling  of  soil  moisture  was  based  on  local-instantaneous  descriptions  of 
hydrologic  processes  such  as  infiltration,  evapotranspiration,  and  recharge.  The  model  describes 
these  inflows  and  outflows  from  the  vadose  zone  as  a  function  of  the  local-instantaneous  soil 
moisture.  These  models  can  be  used  to  describe  the  local  dynamics  of  soil  moisture  or  they  can 
be  integrated  over  spatial  and  temporal  probability  distributions  to  describe  the  large-scale,  long¬ 
term  behavior. 

The  geomorphic  model  used  is  a  fully-distributed  numerical  model.  This  model  was 
extended  to  include  a  more  detailed  hydrologic  description  of  subsurface  flow.  As  a  result,  it  can 
describe  the  occurrence  and  distribution  of  saturated  areas,  their  relation  to  runoff  production, 
and  their  influence  on  topographic  form  and  dynamics. 


4  Main  Results 


4.1  Statistical  Analysis  of  Soil  Moisture  in  Space  and  Time 

The  Illinois  dataset  was  used  to  gain  a  better  understanding  of  the  variability  of  soil  moisture  in 
space  and  time  at  a  relatively  large  scale.  This  analysis  considered  soil  saturation  (the  soil 
moisture  divided  by  the  porosity)  measured  0-10  cm  below  the  ground  surface.  The  observations 
at  each  gage  were  treated  as  independent  and  identically  distributed  random  variables,  which  is 
supported  by  the  low  correlations  observed  between  the  measurements  at  different  gages.  Under 
this  assumption,  the  soil  moisture  can  be  fully  described  by  a  single  probability  density  function 
(PDF).  We  hypothesized  that  this  spatial  PDF  conforms  to  an  Erlang  distribution  and  tested  this 
hypothesis  with  the  Kolmogorov-Smimov  (KS)  test.  This  hypothesis  is  acceptable  for  89%  of 
the  sample  distributions  at  the  10%  level  of  significance.  Furthennore,  we  hypothesized  that  the 
shape  parameter  of  the  Erlang  distribution  remains  constant  for  the  spatial  distributions.  If  this 
hypothesis  holds,  then  the  soil  moisture  patterns  become  more  variable  in  space  when  they  are 
wetter  on  average.  The  KS  test  accepts  the  hypothesis  that  the  samples  are  drawn  from  an  Erlang 
distribution  with  a  constant  shape  parameter  for  94%  of  the  samples  at  the  10%  level.  A  higher 
percentage  of  samples  can  be  accepted  for  this  more  restrictive  case  because  the  KS  test 
measures  only  the  largest  deviation  between  the  sample  distribution  and  the  hypothesized 
distribution,  whereas  the  parameter  estimation  procedure  minimizes  the  error  for  the  distribution 
as  a  whole. 

The  Illinois  data  were  also  used  to  investigate  the  temporal  variability  of  the  spatial  mean 
soil  saturation.  Again,  independence  between  spatial  mean  soil  saturation  at  different  times  was 
assumed.  Two  types  of  analyses  were  performed.  First,  the  distribution  of  mean  soil  saturation 
was  considered  within  a  given  year.  From  this  perspective,  the  parameters  of  the  distribution 
(e.g.,  space-time  mean  soil  saturation)  can  be  specified  for  each  year,  but  the  within-year 
variability  is  described  by  the  distribution.  We  hypothesized  that  the  PDF  conforms  to  a  beta 
distribution.  At  the  10%  level  of  significance,  100%  of  the  temporal  distributions  passed  the  test. 
It  was  also  observed  that  the  standard  deviation  of  the  temporal  distribution  remains 
approximately  constant.  Using  the  KS  test  again,  90%  of  the  samples  were  found  to  conform  to 
a  beta  distribution  with  constant  standard  deviation  at  the  10%  level  of  significance.  The  second 
analysis  considered  the  distribution  for  spatial  mean  soil  saturation  when  all  temporal  variability 


is  included.  This  distribution  includes  both  within-year  and  inter-annual  variability.  Thus,  this 
distribution  describes  the  region’s  climate.  Again,  we  hypothesized  that  a  beta  distribution 
characterizes  the  distribution,  and  this  hypothesis  was  accepted  at  the  10%  level.  Figure  1 
compares  the  sample  temporal  distribution  of  the  spatial  mean  soil  saturation  with  the  beta 
distribution.  A  close  similarity  between  the  two  cumulative  distribution  functions  is  observed. 
The  fact  that  the  beta  distribution  characterizes  both  the  within-year  variability  and  the  total 
variability  is  not  surprising.  The  mean  standard  deviation  for  the  within-year  distribution  is 
0.162  and  the  standard  deviation  for  the  climatic  distribution  is  0.173.  The  similarity  of  these 
two  numbers  implies  that  the  bulk  of  the  variability  of  soil  saturation  in  time  arises  from  within- 
year  variations. 


Figure  1.  Temporal  cumulative  distribution  function  of  spatial  mean  soil  saturation  from  the 
Illinois  dataset  along  with  a  beta  distribution 


Probabilistic  modeling  of  soil  saturation  is  simplified  because  both  the  spatial 
distribution’s  shape  parameter  and  the  temporal  distribution’s  standard  deviation  remain 
relatively  constant  in  time.  In  particular,  if  the  mean  value  of  either  distribution  is  known,  then 
all  characteristics  of  the  PDF  can  be  determined. 

The  correlations  between  the  spatial  mean  soil  moisture  and  the  spatial  mean 
precipitation  and  spatial  mean  potential  evapotranspiration  (PET)  were  also  investigated.  The 
correlation  between  the  mean  soil  saturation  and  precipitation  was  0.16,  which  suggests  that  the 
soil  saturation  tends  to  be  larger  when  the  region  receives  more  precipitation  (as  one  would 


expect).  The  correlation  is  low  because  the  precipitation  rate  is  highly  variable  in  time  whereas 
the  soil  saturation  is  less  variable.  The  correlation  between  the  soil  saturation  and  PET  is  -0.54. 
The  negative  value  implies  that  the  soil  saturation  tends  to  be  low  when  the  PET  is  high.  This 
result  is  also  expected  because  PET  influences  the  loss  of  soil  moisture  through 
evapotranspiration.  Understanding  these  correlations  allows  consideration  of  the  joint 
distributions  of  soil  moisture  and  precipitation  and  soil  moisture  and  PET  in  a  probabilistic 
model.  A  more  detailed  account  of  this  data  analysis  can  be  found  in  Niemann  and  Eltahir 
(2004)  and  Niemann  and  Eltahir  (in  prep.). 


4.2  Modeling  Soil  Moisture  and  Regional  Water  Balance 

The  statistical  analysis  of  soil  moisture  in  the  previous  sections  was  used  to  develop  a  simple 
model  for  soil  moisture  and  regional  water  balance  components.  The  general  approach  of  this 
model  is  to  divide  the  region  into  recharge  and  discharge  locations.  Recharge  locations  can  be 
saturated  or  unsaturated  but  always  have  hydraulic  gradients  oriented  to  allow  recharge  to  the 
groundwater.  Discharge  locations  are  assumed  to  be  saturated  and  have  gradients  oriented  to 
allow  groundwater  discharge  at  the  surface.  The  fluxes  in  and  out  of  the  vadose  zone  are  first 
described  at  the  local  and  instantaneous  scales  for  the  recharge  and  discharge  locations.  Then, 
they  are  integrated  over  the  spatial  and  temporal  distributions  of  soil  saturation,  precipitation, 
and  PET  to  detennine  relations  between  the  space-time  means  of  precipitation,  PET,  soil 
saturation,  evapotranspiration,  surface  runoff,  and  groundwater  discharge. 

The  model  for  infiltration  rate  F  can  be  written  as: 

F  \a(\- s)  +  Kh  Recharge  Locations 
0  Discharge  Locations 


where  a  is  an  infiltrability  parameter,  s  is  the  soil  saturation,  and  Kh  is  the  saturated  hydraulic 
conductivity.  Similar  expressions  have  been  used  previously  (e.g.  Rodriguez-Iturbe  et  ah,  1999), 
but  this  expression  allows  infiltration  for  saturated  soil.  Given  this  infiltration  model,  the  local- 
instantaneous  surface  runoff  R  can  be  written: 


j P-F  if  P>F 
\  0  if  P  <  F 


(2) 


where  P  is  the  local-instantaneous  precipitation  rate.  Coupling  Equations  (1)  and  (2)  shows  that 
runoff  can  occur  because  a  and  Ki,  are  small  (Horton  Runoff)  or  s  is  large  (Dunne  Runoff). 

Two  losses  from  the  vadose  zone  are  allowed.  The  first  is  recharge  to  the  aquifer.  The 
recharge  G  can  be  modeled  as: 

[  Khsy  Recharge  Locations 

G  =  <  (3) 

(  0  Discharge  Locations 


where  /  is  an  exponent  that  depends  on  soil  type.  This  recharge  expression  can  be  derived  from 
Darcy’s  law  if  the  elevation  change  dominates  the  hydraulic  gradient  and  the  Campbell  (1974) 
equation  is  used  to  describe  unsaturated  hydraulic  conductivity.  The  second  loss  is  the 
evapotranspiration  from  the  soil.  If  the  soil  saturation  is  above  the  wilting  point  /?,  it  is  assumed 
that  the  evapotranspiration  is  equal  to  the  PET.  Below  the  wilting  point,  the  evapotranspiration 
depends  on  the  soil  saturation.  In  particular: 


f Eps/J3  ifO  <s<p 
{  Ep  if  /?  <  s  <  1 


(4) 


Similar  approaches  have  been  used  previously  to  describe  the  losses  from  the  vadose  zone  (e.g., 
Dooge  et  ah,  1999),  and  the  combined  loss  function  described  by  Equations  (3)  and  (4)  is 
supported  by  observations  from  Illinois  (Salvucci,  2001). 

After  integrating  over  the  spatial  distributions,  analytical  expressions  were  obtained  for 
the  spatial  mean  surface  runoff,  evapotranspiration,  and  groundwater  recharge  as  a  function  of 
the  spatial  mean  soil  saturation.  Then,  by  integrating  over  the  temporal  distributions  of  the 
spatial  means,  similar  relationships  were  detennined  for  the  space-time  mean  values.  By 
applying  equilibrium  conditions  for  the  groundwater  and  the  system  as  a  whole,  the  space-time 
mean  soil  saturation  can  be  determined  from  the  space-time  mean  precipitation,  PET,  and 
saturated  hydraulic  conductivity.  These  equilibrium  conditions  hold  in  approximation  at  the 
annual  and  climatic  time  scales.  Figure  2a  compares  the  annual  runoff  determined  from  the 
model  with  annual  runoff  observations  from  Illinois.  The  figure  shows  that  the  model  captures 
the  basic  tendency  of  the  observations.  In  this  model,  only  precipitation  was  varied  as  an  input, 
so  the  scatter  observed  in  this  plot  cannot  be  reproduced.  The  figure  also  shows  the  relative 
contributions  of  surface  runoff  and  groundwater  discharge  to  the  total  runoff  (distinguished  by  a 
dotted  line).  Figure  2b  graphs  the  dependence  of  the  runoff  coefficient  on  the  humidity  index 
De,  which  is  defined  as  the  mean  precipitation  divided  by  the  mean  PET.  Several  empirical 


expressions  have  been  proposed  that  describe  the  observed  variation  of  the  runoff  coefficient 
between  different  regions  or  years  as  a  function  of  the  humidity  index  (Dooge,  1992).  The 
runoff  coefficient  calculated  by  the  model  actually  depends  on  two  humidity  indices:  DE  and  DK, 
where  DK  is  the  mean  precipitation  divided  by  the  saturated  hydraulic  conductivity.  The  blue 
lines  in  Figure  2b  demonstrate  that  the  model’s  dependence  on  DE  is  similar  to  Schreiber’s 
empirical  equation.  As  the  humidity  index  becomes  large,  the  evapotranspiration  is  not  restricted 
by  the  availability  of  water  at  as  many  locations  and  times.  As  a  result,  the  evapotranspiration 
rate  approaches  the  PET  rate,  which  constrains  its  ability  to  transport  water.  This  restriction 
leads  to  an  increase  in  mean  soil  saturation,  which  promotes  Dunne  runoff  and  thus  higher  runoff 
coefficients. 


Figure  2.  (a)  Comparison  of  observed  annual  runoff  from  Illinois  with  model  predictions,  and 
(b)  comparison  of  modeled  runoff  coefficient  with  three  empirical  equations 

Another  key  result  from  the  model  is  the  nature  of  the  relationship  between  the  soil 
saturation  and  the  hydrologic  and  physical  characteristics  of  the  region.  Figure  3  shows  the 
derived  relationship  between  the  mean  soil  saturation  and  the  two  humidity  indices  DE  and  DE. 
As  the  precipitation  becomes  larger  with  respect  to  either  PET  or  saturated  hydraulic 
conductivity,  the  mean  soil  saturation  increases  as  expected.  The  nature  of  this  increase  can  be 
classified  into  three  regimes.  When  DE  is  low,  the  value  of  the  hydraulic  conductivity  has  very 
little  influence  on  the  mean  soil  saturation,  so  the  region  is  DE  controlled.  In  this  regime, 
evapotranspiration  is  the  dominant  loss  from  the  vadose  zone.  Similarly,  when  DE  is  small,  the 


hydraulic  conductivity  is  large,  so  groundwater  recharge  becomes  a  major  loss.  However,  even 
for  very  low  Dk  values,  the  influence  of  evapotranspiration  is  observed.  In  the  third  regime, 
evapotranspiration  is  the  dominant  loss,  but  groundwater  recharge  is  also  significant.  This  zone 
is  considered  DE  dominated.  See  Niemann  and  Eltahir  (2004)  for  more  details  on  this  subject. 


Space-Time  Mean  Soil  Saturation 


Figure  3.  Estimation  of  space-time  mean  soil  saturation  from  precipitation,  PET,  and  hydraulic 
conductivity 

The  modeling  approach  was  also  used  to  understand  the  impact  of  a  region’s  physical  and 
statistical  characteristics  on  the  mean  soil  saturation  and  water  balance  components.  In 
particular,  we  examined  the  sensitivity  of  soil  moisture,  surface  runoff,  and  groundwater 
discharge  to  changes  in  precipitation  and  PET.  The  sensitivity  is  defined  as  the  percentage 
change  in  a  variable  when  the  mean  precipitation  or  PET  is  changed  by  one  percent  (Schaake, 
1990;  Sankarasubramanian  et  ah,  2001).  Analysis  of  the  sensitivities  is  useful  not  only  for 
characterizing  soil  moisture  but  also  for  forecasting  the  vulnerability  of  a  region’s  water  cycle  to 
climate  changes.  Figure  4  plots  the  sensitivity  of  mean  soil  saturation,  surface  runoff,  and 
groundwater  discharge  to  changes  in  PET.  The  sensitivities  are  calculated  for  regions  that  differ 
only  in  their  wilting  point  (f:J)  or  their  temporal  standard  deviation  of  soil  saturation  (cr).  The 
figure  shows  that  regions  with  higher  wilting  points  have  less  sensitive  mean  soil  saturations. 
When  the  wilting  point  is  large,  the  evapotranspiration  is  more  limited  by  moisture  availability, 
so  the  fluxes  and  soil  moisture  are  less  sensitive  to  changes  in  PET.  The  figure  also  shows  that 
regions  with  more  temporal  variability  of  soil  saturation  have  less  sensitive  surface  and 
groundwater  runoff  and  slightly  more  sensitive  soil  saturation.  Temporal  variability  of  soil 


saturation  smoothes  out  thresholds  in  the  partitioning  of  local-instantaneous  precipitation.  Thus, 
more  variability  in  soil  moisture  results  in  smoother  relations  and  less  sensitivity  of  the  fluxes.  A 
more  detailed  description  of  the  sensitivity  of  soil  moisture  can  be  found  in  Niemann  and  Eltahir 
(in  prep.). 
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Figure  4.  Sensitivity  of  space-time  mean  soil  saturation  and  hydrologic  fluxes  to  changes  in 
PET  as  the  wilting  point  (fJ)  and  the  temporal  standard  deviation  of  soil  saturation  (a)  are 
changed. 

4.3  Dependence  of  Soil  Moisture  on  Regional  Characteristics 

The  previous  subsections  described  some  statistical  properties  of  soil  moisture  and  presented  a 
model  that  relates  soil  moisture  to  precipitation,  PET,  and  losses  from  the  vadose  zone.  In  this 
subsection,  the  observed  dependence  of  soil  saturation  on  other  physical  characteristics  of  a 
region  are  analyzed. 

Soil  saturation  patterns  are  thought  to  have  a  complex  dependence  on  regional 
topography,  soil  type,  and  vegetation.  The  SGP  97  soil  saturation  data  have  a  high  degree  of 
variability  in  both  space  and  time.  In  this  project,  we  hypothesized  that  this  variability  can  be 
well  described  by  a  small  number  of  underlying  independent  spatial  patterns.  These  patterns  are 
thought  to  reflect  certain  characteristics  of  the  region  such  as  percent  sand,  percent  clay,  surface 
roughness,  bulk  density,  and  vegetation  water  content.  However,  the  connections  between  these 
patterns  and  the  regional  characteristics  may  depend  upon  the  spatial  and  temporal  scales 
considered. 


Empirical  Orthogonal  Function  (EOF)  and  Principal  Component  (PC)  analysis  was  used 
to  describe  the  total  variability  from  a  16  day  time  series  of  observations.  The  type  of  variability 
explained  by  EOF  analysis  depends  on  the  preparation  of  the  data.  We  have  analyzed  the  SGP 
97  in  two  ways.  First,  the  average  soil  saturation  value  was  found  for  the  time  series  at  each  grid 
point  and  removed  to  calculate  temporal  soil  saturation  anomalies.  The  resulting  EOFs  are 
spatial  patterns  (maps),  explaining  the  variability  the  grid  points  relative  to  their  long-term 
averages.  Second,  the  data  was  analyzed  in  spatial  anomaly  form  by  removing  the  spatial 
average  on  each  day  from  every  grid  point.  When  the  data  is  placed  into  spatial  anomaly  form, 
the  EOFs  are  maps  explaining  the  variability  of  each  grid  point  relative  to  the  spatial  average. 

EOF  analysis  is  based  on  the  notion  that  a  data  matrix  X  (e.g.  spatial  soil  saturation 
anomalies)  can  be  decomposed  into  PCs,  eigenvalues,  and  EOFs.  In  particular: 

X  =  UEVt  (5) 

where  the  columns  of  U  contain  the  PCs,  the  diagonal  values  of  E  contain  the  squares  of  the 
eigenvalues,  and  the  rows  of  V  contain  the  EOFs.  The  principle  components  relate  to  the 
temporal  behavior,  the  eigenvalues  are  used  to  rank  the  significance  of  the  EOF/PC  pairs,  and  the 
EOFs  relate  to  the  spatial  patterns.  Soil  saturation  anomalies  were  analyzed  to  find  EOFs  by  first 
multiplying  the  anomalies  matrix  by  its  transpose  to  find  the  autocorrelation  matrix.  This 
autocorrelation  matrix  can  be  analyzed  to  find  eigenvectors  (corresponding  to  the  PCs)  and 
eigenvalues.  To  remove  any  nonphysical  amplitude  that  results  from  the  eigenanalysis,  the  PCs 
were  standardized  to  have  a  mean  of  zero  and  a  standard  deviation  of  one.  By  multiplying  the 
standardized  PCs  with  the  original  anomaly  matrix,  one  can  find  the  EOFs,  which  have  the  same 
units  as  the  original  data.  The  variance  explained  by  each  of  the  EOFs  is  equal  to  the 
corresponding  eigenvalue  divided  by  the  sum  of  all  the  eigenvalues.  EOF/PC  pairs  are 
considered  significant  if  there  is  no  overlap  of  their  error  bar  plots.  Error  bars  are  detennined  by 
the  fraction  of  variance  explained  (A)  plus  or  minus  (A(2 /N)  )/2  where  N  is  equal  to  the  number 
of  independent  samples  (taken  as  the  number  of  days  in  the  time  series). 

When  the  soil  saturation  data  (in  temporal  anomaly  form)  were  analyzed,  five  significant 
EOF/PC  pairs  were  found.  Together  they  explain  94%  of  the  variance  in  the  original  data,  and 
50%  of  the  variance  is  explained  by  temporal  EOF  1  alone.  When  the  data  were  analyzed  a 
second  time  in  spatial  anomaly  form,  five  significant  EOF/PC  pairs  were  also  found  which 
together  explain  94%  of  the  variance.  However,  the  first  spatial  EOF  explains  only  61%  of  the 


variance.  Spatial  EOF  1  is  shown  in  Figure  5  along  with  the  soil  saturation  pattern  from  June  18 
and  the  pattern  of  percent  sand. 


Soil  Saturation  6/18  Spatial  EOF  1  %  Sand 


Figure  5:  A  sample  map  of  soil  saturation  data,  spatial  EOF  1,  and  percent  sand  for  visual 
comparison. 

EOF  analysis  can  be  performed  in  a  variety  of  ways,  depending  on  one’s  objectives.  We 
have  found  the  EOFs  and  PCs  from  a  data  matrix  that  contains  only  the  soil  saturation  fields. 
Kim  and  Barros  (2002)  also  generated  EOFs  and  PCs  using  the  SGP  97  data.  Their  EOF  analysis 
was  perfonned  using  a  matrix  with  soil  saturation,  vegetation  water  content,  percent  sand  and 
topography.  Their  EOFs  illustrate  patterns  that  captured  the  interdependencies  of  each  of  these 
variables.  The  EOFs  we  calculated  describe  the  spatial  and  temporal  variability  of  soil  saturation 
in  terms  of  a  small  number  of  independent  spatial  patterns. 

The  correlation  coefficient  is  a  common  way  to  determine  whether  two  variables  have  a 
linear  dependency.  We  have  used  it  to  measure  how  closely  our  EOFs  are  related  to  the  physical 
characteristics  of  the  region.  When  the  analysis  was  performed  on  spatial  EOF  1 ,  it  was  found  to 
have  a  correlation  coefficient  of  -0.32  with  percent  sand,  0.04  with  percent  clay,  -0.07  with 
surface  roughness,  -0.20  with  vegetation  water  content,  and  -0.20  with  bulk  density.  The  high 
(negative)  correlation  with  percent  sand  is  an  expected  result.  For  example,  visual  inspection 
shows  that  the  southern  portion  of  the  region  tends  to  have  high  percent  sand  values  and  low 
EOF  1  values  (Figure  5). 


The  SGP  97  data  itself  was  correlated  on  a  daily  basis  with  regional  characteristics  by 
Kim  and  Barros  (2002).  They  found  that  on  ‘wet’  days,  soil  saturation  was  most  highly 
correlated  with  percent  sand,  while  on  ‘dry’  days  it  was  most  highly  correlated  with  percent  clay. 
They  also  found  that  each  of  the  ancillary  data  correlations  increased  with  aggregation.  In 
contrast,  our  research  has  found  a  set  of  EOF  maps  that  efficiently  describe  the  variability  of  soil 
saturation  for  the  duration  of  SGP  97  experiment.  The  most  significant  of  these  maps  is  highly 
correlated  with  percent  sand,  but  not  well  correlated  with  percent  clay. 

The  correlation  analysis  described  above  was  performed  with  the  data  at  the  original 
pixel  size  of  0.64  km  .  Thus,  the  analysis  includes  spatial  variability  from  a  wide  range  of  spatial 
scales  (from  0.64  km"  to  nearly  the  size  of  the  domain).  To  detennine  which  variables  are 
related  to  the  spatial  pattern  of  EOF  1  at  different  scales,  we  repeated  the  correlation  analysis  on 
aggregated  data.  Both  the  EOFs  and  the  regional  characteristic  data  were  aggregated  from  an 
original  cell  size  of  0.64km  up  to  a  cell  size  of  1024km  .  Figure  6  shows  the  results  of  this 
analysis.  The  correlation  of  EOF  1  with  percent  sand,  percent  bulk  density,  and  vegetation  water 
content  increases  slightly  with  increasing  aggregation.  The  correlation  of  EOF  1  with  percent 
clay  starts  out  slightly  positive  but  turns  slightly  negative  as  the  aggregation  scale  increases.  A 
dramatic  increase  in  the  correlation  of  EOF  1  and  surface  roughness  is  observed  with  increasing 
scale.  A  visual  inspection  on  aggregated  maps  of  EOF  1  and  surface  roughness  revealed  that 
locations  with  low  surface  roughness  also  have  low  values  for  EOF  1.  The  observed  positive 
correlation  at  large  scales  may  be  associated  with  the  fact  that  smooth  topography  increases  the 
potential  for  evapotranspiration.  Smooth  topography  is  more  uniformly  oriented  towards  the  sun 
and  therefore  receives  more  radiation. 


Spatial  E0F1  Correlations 


Figure  6.  Correlation  coefficients  of  spatial  EOF  1  with  regional  characteristics  at  a  range  of 
scales. 


Another  project  objective  is  to  understand  the  relationship  between  the  scaling  properties 
of  topography  and  those  of  soil  saturation.  Numerous  indices  are  available  to  estimate  the 
relative  soil  wetness  of  a  location  from  topographic  information.  These  indices  usually 
determine  the  soil  saturation  (or  a  proxy)  as  a  deterministic  function  of  topography  and  other 
data  (e.g.,  O’Loughlin,  1986).  One  such  index  is  the  so-called  Topindex  (Beven  and  Kirkby, 
1979),  which  we  used  to  analyze  the  relationship  between  the  scaling  properties  of  topography 
and  those  of  soil  saturation.  A  well-known  implication  of  the  scaling  invariance  of  topography  is 
the  so-called  slope-area  relationship  (Veneziano  and  Niemann,  2000),  which  states  that  the  slope 
of  a  channel  decreases  as  the  topographic  area  that  drains  through  the  channel  increases.  In 
particular,  S  «  bA  "  where  S  is  channel  slope,  A  is  the  drained  area,  b  is  a  coefficient,  and  6  is  an 
exponent  that  typically  falls  within  0.2  <0  <0.6  (Veneziano  and  Niemann,  2000).  Figure  7 
plots  the  Topindex  as  a  function  of  drained  area  (“sub-basin  area”).  The  figure  shows  three 
distinct  segments  in  this  relationship.  When  the  area  draining  through  the  point  is  small,  the 
point  is  located  on  a  hillslope  and  the  Topindex  exhibits  little  dependence  on  the  drained  area. 
As  the  drained  area  increases,  the  points  are  located  in  unchannelized  valleys  and  in  ephemeral 
streams  and  a  power  law  is  observed.  At  very  large  drained  areas,  the  points  fall  within 
perennial  channels  and  the  Topindex  reaches  a  maximum  saturated  value.  Analytical 


expressions  were  also  derived  that  describe  the  shape  of  this  scaling  relationship  in  terms  of  the 
parameters  of  a  common  landscape  evolution  model  (from  Moglen  and  Bras,  1995). 


Figure  7.  Scaling  of  Topindex  with  sub-basin  or  drained  area 


4.4  Topographic  Scaling  Invariance  and  Its  Origin 

Another  line  of  research  was  an  investigation  of  the  origin  and  nature  of  topographic  scaling 
invariance.  Understanding  topographic  scaling  is  important  if  this  condition  is  to  be  used  to  infer 
fine-scale  soil  moisture  in  part  from  coarse-scale  topographic  information.  As  part  of  this  line  of 
research,  the  scaling  properties  of  a  series  of  laboratory-generated  topographies  were  analyzed 
(Hasbargen  and  Paola,  2000).  These  topographies  were  generated  under  controlled  rainfall  rates, 
uplift  rates,  and  substrate  characteristics.  Under  most  cases  considered,  the  substrate  was  kept 
saturated,  so  the  runoff  production  was  constant  in  space  and  time.  One  of  the  key  results  from 
these  laboratory  experiments  was  an  observed  mobility  of  the  channels  even  when  the 
denudation  rate  balanced  the  uplift  rate  on  average.  Under  such  conditions,  most  numerical 
landscape  evolution  models  produce  temporally  invariant  surfaces  (Hasbargen  and  Paola,  2000). 
We  examined  the  scaling-invariance  of  topographies  generated  under  these  conditions  and 
compared  it  with  the  scaling-invariance  of  typical  natural  basins.  Our  analysis  identified  a 
quantitative  scaling  property  that  coincided  with  the  mobility  of  channels.  In  particular,  small 
tributaries  near  large  channels  typically  have  lower  slopes  than  expected  from  the  mean  slope- 
area  relationship  for  the  experiments  that  exhibited  mobility.  Figure  8  demonstrates  this 
tendency.  The  left  column  of  plots  shows  the  slope-area  coefficient  as  function  of  the  ratio  of  a 
nearby  channel’s  drained  area  and  the  local  drained  area.  For  the  Madulce  basin  in  California 


the  coefficient  is  independent  of  this  ratio  of  areas.  For  the  laboratory  experiments  (both  5  and 
6),  the  coefficient  decreases  slightly  as  the  ratio  increases.  This  decrease  implies  that  the 
presence  of  nearby  large  channels  results  in  smaller  slopes  for  small  channels.  The  right  column 
of  plots  displays  the  same  analysis  but  the  drained  areas  are  determined  after  moving  specified 
distances  upstream.  The  dependence  disappears  with  increased  distance  in  Experiment  5,  which 
had  less  lateral  movement,  but  it  remains  visible  for  Experiment  6,  which  had  more  lateral 
movement.  This  property  and  other  scaling  characteristics  of  laboratory  and  natural  basins  are 
described  in  Niemann  and  Hasbargen  (in  review). 
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Figure  8.  Dependence  of  secondary  tributary  slopes  on  primary  tributary  slopes. 


4.5  Relationship  between  Subsurface  Flow  and  Topographic  Dynamics 

The  final  line  of  research  examined  the  impact  of  subsurface  flow  on  the  evolution  of 
topographic  scaling  and  the  associated  evolution  of  hydrologic  processes.  The  general  approach 
used  to  study  these  issues  was  the  development  of  an  improved  numerical  landscape  evolution 
model.  The  main  improvement  lies  in  the  description  of  the  hydrologic  processes.  In  particular, 
infiltration  and  groundwater  dynamics  have  been  included  when  determining  runoff  production 
and  streamflow. 

The  erosion  model  used  is  a  variant  of  the  GOLEM  landscape  evolution  model  (Tucker 
and  Slingerland,  1997).  The  governing  equation  under  the  detachment-limited  condition  can  be 
written: 

—  =  U  —  kQmSn  +  DV2z  (6) 

8t 

where  z  is  elevation,  t  is  time,  U  is  uplift  rate,  k  is  an  erodability  coefficient,  Q  is  water 
discharge,  S  is  channel  slope,  m  and  n  are  exponents,  and  D  is  the  hillslope  diffusivity.  In  most 
landscape  evolution  models,  the  discharge  rate  Q  is  calculated  from: 

Q  =  PA  (7) 


where  Pe  is  an  effective  precipitation  rate  and  A  is  the  drained  area.  Pe  is  spatially  uniform, 
which  disallows  Dunne  runoff  production  and  groundwater  discharge. 

In  the  improved  model,  runoff  production  and  groundwater  discharge  are  detennined 
dynamically.  Precipitation  is  simulated  as  a  series  of  discrete  storm  events  using  the  rectangular 
pulse  model  (Eagleson,  1978;  Tucker  and  Bras,  2000),  which  includes  three  independent  random 
variables:  rainfall  intensity  P,  rainfall  duration  7)>,  and  inter-storm  period  To-  Runoff  and 
infiltration  are  determined  using  a  model  similar  to  the  one  described  earlier  in  this  report. 
However,  in  these  preliminary  simulations,  infiltration  and  recharge  rates  are  assumed  to  be 
independent  of  the  soil  saturation  for  simplicity.  The  groundwater  is  simulated  using  the  Dupuit 
approximation,  which  assumes  that  the  hydraulic  gradient  coincides  with  the  water  table  and  the 
groundwater  flow  is  mainly  horizontal.  Under  these  assumptions  the  water  table  level  h  is 
governed  by  the  expression: 


~d\h2)  d2(h2) 
8x2  +  8y2 


+  F 


(8) 


where  x  and  y  are  the  horizontal  coordinates  and  Sy  is  the  specific  yield.  This  expression 
assumes  that  the  underlying  impermeable  layer  is  horizontal  (for  simplicity  of  notation),  but  the 
numerical  model  can  include  spatial  variations  in  the  elevation  of  this  layer  as  well. 
Groundwater  discharge  occurs  at  locations  where  the  water  table  meets  the  land  surface.  The 
total  streamflow  is  the  sum  of  the  surface  runoff  and  groundwater  discharge. 

Figure  9  shows  a  series  of  simulated  precipitation  events  along  with  the  streamflow 
produced  by  the  model.  The  streamflow  exhibits  spikes,  which  are  associated  with  the  surface 
runoff  from  precipitation  events,  and  base  flow,  which  is  associated  with  groundwater  discharge. 


Figure  9.  Simulated  precipitation  and  the  associated  streamflow  at  the  simulated  region’s  outlet 


The  model  has  been  used  to  investigate  the  impact  of  hydrologic  processes  on  bedrock 
detachment,  sediment  transport,  topographic  form,  and  the  evolution  of  the  hydrologic  processes 
themselves.  Figure  10  shows  the  proportion  of  a  simulated  topography  that  is  saturated  for 
events  with  different  durations  and  return  periods.  The  two  lines  show  cases  where  the 
topography  has  evolved  with  and  without  substantial  groundwater  flow  (the  groundwater  flow  is 
controlled  by  changing  the  saturated  hydraulic  conductivity).  Interestingly,  the  topography  with 
substantial  groundwater  flow  evolves  toward  a  form  that  promotes  saturated  areas  particularly 
for  longer  duration  storms.  The  same  topography  also  evolves  toward  a  form  that  promotes 


runoff,  especially  from  less  frequent  events  (i.e.  large  return  periods).  These  results  hint  at  a 
potential  feedback  loop  whereby  patterns  of  runoff  production  and  groundwater  discharge 
modify  the  topographic  surface  to  enhance  their  efficiency  for  certain  types  of  events. 


Rainfall  Duration  (hours)  Return  Period  (years) 


Figure  10.  Impact  of  groundwater  runoff  on  the  portion  of  a  basin  that  saturates  during  events 
with  different  durations  (a)  and  return  periods  (b) 

5  Impacts 

A  wide  variety  of  Army  activities  could  benefit  from  soil  moisture  information  at  fine  spatial  and 
temporal  resolutions.  For  example,  soil  moisture  information  can  be  used  to  assess  the  mobility 
of  troops  and  vehicles  during  combat,  including  travel  time  calculations  and  optimal  travel  path 
determinations.  Soil  moisture  information  can  also  be  used  to  estimate  vegetation  stem  spacing, 
which  impacts  mobility.  Another  activity  that  depends  on  soil  moisture  is  land-mine  detection. 
This  application  would  require  soil  moisture  information  at  very  fine  resolutions.  Land 
management  practices  can  also  benefit  from  a  better  understanding  of  soil  moisture  patterns  and 
dynamics.  For  example,  the  location  and  timing  of  training  exercises  could  be  selected  in  part  to 
avoid  excessive  degradation  of  training  lands. 

The  research  conducted  in  this  project  provides  an  improved  base  of  knowledge  that  will 
be  used  to  develop  soil  moisture  interpolation  procedures  that  include  the  role  of  topography. 
Topographic  infonnation  is  available  for  nearly  all  regions  of  the  globe  at  a  resolution  that  is 
consistent  with  many  Army  activities.  If  topographic  information  can  be  quantitatively  related  to 
soil  moisture  patterns  in  an  efficient  manner,  this  information  can  be  used  to  estimate  fine-scale 
soil  moisture  patterns  from  available  data.  In  addition,  a  better  understanding  of  the  relationships 


between  soil  moisture  and  climatic  variables  such  as  precipitation  and  PET  provides  a  better 
ability  to  forecast  changes  in  soil  moisture  with  time. 

The  project  will  also  benefit  a  range  of  scientific  fields.  A  new  model  for  regional  water 
balance  was  presented  that  can  be  used  to  assess  the  impact  of  climatic  fluctuations  and  changes 
on  the  water  cycle  of  a  region.  In  addition,  the  new  process-based  landscape  evolution  model  is 
expected  to  have  impacts  on  various  geomorphic  applications.  For  example,  the  model  could 
eventually  be  used  to  assess  the  fate  of  pollutants  on  Army  lands  over  long  periods  of  time  (when 
topography  must  be  considered  dynamic).  Finally,  the  research  conducted  in  this  project  is 
expected  to  benefit  the  estimation  and  simulation  of  soil  moisture  for  civilian  hydrologic 
applications  such  as  water  resources  planning  and  irrigation/agricultural  policy  development. 
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